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ABSTRACT 

We give an algorithm for the economical calculation of angles and actions for stars 
in axisymmetric potentials. We test the algorithm by integrating orbits in a realistic 
model of the Galactic potential, and find that, even for orbits characteristic of thick- 
disc stars, the errors in the actions are typically smaller than 2 percent. We describe a 
scheme for obtaining actions by interpolation on tabulated values that significantly ac- 
celerates the process of calculating observables quantities, such as density and velocity 
moments, from a distribution function. 



1 INTRODUCTION 



When electronic computers first became widely available, it 
was discovered that orbits in typical axisymmetric galactic 
potentials usually admit three isolating integrals of motion 
(Henon & Heiles 1964; Ollongren 1965). Consequently, by 
Jeans' theorem, the distribution functions (dfs) of equilib- 
rium axisymmetric galaxies should be functions of three in- 
tegrals of motion. Unfortunately, analytic forms of all three 
integrals are known only for exceptional potentials, so the 
few three-dimensional galaxy models in the literature that 
have a known df (e.g. Rowley 1988) employ only the classi- 
cal energy and angular-momentum integrals E and L z , and 
therefore lack generality. 

The action integrals J r , J z and L z are particularly use- 
ful constants of motion (e.g. Binney 2012), and we have 
previously argued the merits of models in which the dis- 
tribution function is an analytic function of J r , J z and L z . 
To take advantage of these models one should be able to 
evaluate economically the actions of a star from its conven- 
tional phase-space coordinates (x, v). To date we have used 
two techniques for evaluating actions: (i) torus construction 
(Kaasalainen & Binney 1994; Binney & McMillan 2011) and 
(ii) the adiabatic approximation (Binney 2010; Binney & 
McMillan 2011; Schonrich & Binney 2012). Torus construc- 
tion is a general and rigorous technique and for some appli- 
cations it is the technique of choice (e.g. McMillan & Binney 
2012). For other applications it is inconvenient because it 
delivers (x, v) as functions of the actions and angles, rather 
than the actions and angles as functions of (x,v). 

The adiabatic approximation delivers actions and an- 
gles as functions of (x, v) but it is reasonably accurate only 
for stars that stay close to the Galaxy's mid-plane. Here 
we introduce a different approximate way to obtain actions, 
which, though still approximate, is more accurate than the 
adiabatic approximation and is valid for stars that move far 
from the mid plane. 



2 THE ALGORITHM 

Our algorithm is based on the idea that the Galaxy's gravi- 
tational potential is similar to a Stackel potential - for a de- 
tailed description of the latter see de Zeeuw (1985). Stackel 
potentials for oblate bodies are framed in terms of prolate 
confocal coordinates. The latter are defined by the distance 
2A between the foci of the coordinate curves. These foci lie 
at R = and z — ±A, where (R, z, <j>) is a system of cylindri- 
cal polar coordinates. Following Binney & Tremaine (2008; 
hereafter BT08) §3.5.3 we define new coordinates (u, v) by 



R = A sinh u sin v 



z = A cosh u cos -u. 



(1) 



The generating function of the canonical transformation be- 
tween these systems of coordinates is 



S(pn,pz,u,v) = p R R{u,v) + p z z(u,v) 

so from p u = dS/du we have 

p u = A(pn cosh u sin v + p z sinh u cos v) 
p v = A(pR sinh u cos v — p z cosh u sin v) 



(2) 



(3) 



In these coordinates a Stackel potential can be written in 
terms of two functions of one variable, U(u) and V(v), being 
given by 



U(u) - V(v) 
sinh 2 u + sin 2 v 



(4) 



This being so, the Hamilton-Jacobi equation yields (BT08 
eq. 3.249) 



^= J Bsinh 2 u-J 3 -[/H 

2 

JbL [ =Esm 2 v + h + V(v)- 



2A 2 sinh 2 u 



(5) 



where E is the orbit's energy and 73 is a constant of sep- 
aration. These equations make p u (u) and p v (v) functions 
of only their conjugate coordinate, so we can evaluate the 
actions as 
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IT 



dup u (u) ; J z = - 

IT 



r/2 



dwp„(v), (6) 



where tt m in < u max are the roots of p u (u) = and v m i n 
is the root of p v (v) = 0. Note that an orbit's actions are 
independent of any system of coordinates and the subscripts 
r and z on the actions merely remind us that, in a general 
way, J r quantifies oscillations inwards and outwards, while 
J z quantifies oscillations around the equatorial plane. 

In as much as our potential $ is similar to a Stackel 
potential, we have 



(sinh 2 u + sin 2 v)<J>(m, v) ~ U (u) — V(v). 



(7) 



Consequently, we have 

SU = (sinh 2 u + sin 2 v)&(u, v) — (sinh 2 w + sin 2 v)&(u , v) 

~ U{u) - U(u Q ) 
SV = cosh 2 w$(w, 7r/2) — (sinh 2 w + sin 2 u)$(u, «) (8) 

~V(v) - V(tt/2). 

Here Mo is a reference value of u, the choice of which will 
be discussed below, and the right side of the first equation 
appears to be a function of v but its dependence on v will be 
weak unless $ is very unlike a Stackel potential. Similarly, we 
assume that the dependence of the right side of the second 
equation on u is at most weak. Then, given a point (x, v) 
on the orbit we can calculate two constants of motion: 



h + U(u ) ~ h + U(u) - 6U(u) 



E sinh u - 



pi 



Li 



2A 2 2A 2 sinh 2 u 
h + V(n/2) ~ J 3 + V» - <5V» 

__ Ssm u+ _ 



2A 2 sin 2 v 

Now we can evaluate p u for any given w from 



SU(u) 



5V{v). 



(9) 



-,(10) 



J* * B sinh 2 u - [/a + f/Cno) + - , 

so we can evaluate the integral for J r . The integral for J z is 
evaluated in the same way. 

In principle uo can be taken to be any quantity that is 
constant along an orbit, but the accuracy of our work will 
depend on our choosing a value such that the term in the 
definition (8) of SU that contains uo almost completely elim- 
inates the v dependence of the first term in this equation. 
In fact, the natural choice for uo is the location u of the 
minimum with respect to u of SU at fixed v. This minimum 
can be determined before we have specified uo because the 
derivative with respect to u of the first of equations (8) is 
manifestly independent of uo- Physically u is the radial co- 
ordinate of the shell orbit J r = of given values of E and 
Lz. 



2.1 Angle variables 

Equations (3) for the momenta are obtained by solving 
the Hamilton-Jacobi equation for the generating function 
S{u, v, 4>, J r ,Jz,L z ) of the canonical transformation between 
the (u, v, 4>,p u , ■ ■ ■) and the (6 r ,9z,0^, J r , . . .) systems of 
canonical coordinates with S of the form 



Given that S takes this form, we may write 
= J Aup u + J i 



dvp v + 4>L Z 



(12) 



Hence 



9 = — = [ du — + 

8Jr / dJ r 



dv 



dp v 

dJ r 

dp v 



(13) 



6 * = Wz = j dU Wz + j AV dJ; 

dS f dp u f dp v 

9 * = dL M =J du dTz + J dv oTz +4> - 

We obtain the derivatives of p u and p v from the chain rule. 
For example 

dp u _ dpn_dE_ dp u dl 3 
dJ r dE dJ r dh dJ r 

~ dE " dl 3 dJr ' 

where 0, r = dE/dJ r is the radial frequency, so 



(14) 



v/2A 



, sinh 2 u 
du h 



dh 

8J r 



Pu 

du 

Pu 



dv ■ 



pv 



du 



(15) 



A detail possibly worth noting is that we always take p u of 
p v to be given by the positive square root and when con- 
sidering a point in phase space at which p u < we obtain 
the indefinite integrals over u as twice the corresponding in- 
tegral from tt m i n to u max minus the integral from u m i n to 
u with p u taken to be positive. When this procedure is fol- 
lowed for all integrals, the angle variables increase along an 
orbit continuously as they should. 

The derivatives with respect to J r in equation (14) can 
be obtained by observing that by the chain rule the matrix 

(fl r Q. z fltp \ 

dl 3 /dj r dh/dJz dh/dLz (16) 
1 



is the inverse of the matrix 1 

fdJr/dE dJ r /dh OJr/dLz 
dJz/dE dJz/dh OJz/dLz 
V 1 



(17) 



S = S U {U, J r , Jz, L z ) + S V (V, J r , Jz,L z ) + (j>L z 



(11) 



The latter is readily obtained by differentiating equations (6) 
and leads to the definite integrals mentioned in the previous 
paragraph. 

2.2 Interpolation 

To recover the observable properties of a model stellar sys- 
tem at a given spatial point, such as its density p and velocity 
dispersion tensor a 2 , , one has to integrate the distribution 
function over all velocities. These integrals entail large num- 
bers of evaluations of the df, and it is important to keep 

1 Care must be taken with derivatives with respect to L z regard- 
ing whether they arc at constant (E, I3) or (J r , J z ). 
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down the cost of each evaluation. This goal motivates us to 
tabulate the values of J r and J z as functions of the classical 
integrals E, L z and -Z3 + U(uo) or 73 + V(n/2). However, 
I3 + U(uo) proves ill-suited to this task because its numeri- 
cal value varies rapidly as one moves through action space. 
A more convenient constant of motion is 

Er S <fe + it (^dk - ^mW) + SU(U) 

— i?(sinh 2 u — sinh 2 uq). (18) 

At u — uo, which we have chosen to be the minimum of the 
potential that governs the motion in u, E r = p 2 /2A 2 so we 
can think of E r as the energy invested in radial oscillations. 
Consequently, for any values of E and L z , E r vanishes for 
J r = and takes its largest value for J z = and we can 
readily obtain J r and J z by interpolating between the values 
taken by J r and J z at a grid of values of E r . 

In detail we structure the grid in (L z , E, E r ) space as 
follows. The grid points in L z are defined by the angular 
momenta of circular orbits with radii uniformly distributed 
between minimum and maximum radii. For each value of L z 
we adopt as grid points in E the energies 

E t = E C (L Z )+ (5)7 w) 2 , (19) 

where E C (L Z ) is the energy of the circular orbit with angular 
momentum L z and ^v^ ax is slightly smaller than the differ- 
ence between the energy of that orbit and the escape energy 
from its circle. For each such energy we identify uo = u, the 
minimum with respect to u of 

E sinh 2 u - 5U -L 2 Z / (2 A 2 sinh 2 u). (20) 

Then we find the speed v that the star has at this spa- 
tial point and determine the values taken by E r , I3 + 
V(n/2), J r and J z at the phase-space point (x, v) = 
(A sinh(tto), 0, v cos ip, v sin ip) for values of tp uniformly dis- 
tributed in (0,7r/2). With this scheme interpolation errors 
can be kept below ~ 1% with a grid of size 60 x 50 x 50, 
which takes ~ 30 sec to compute on a laptop. 

The present algorithm lends itself to tabulation better 
than the adiabatic approximation because with the present 
algorithm it is straightforward to resort to the algorithm 
whenever actions are required for values of the integrals 
that lie outside the grid. By contrast, when the adiabatic 
approximation is used, values of E z are required for given 
J z and these are hard to obtain beyond the limits of the 
pre-computed table of values of J z for given E z . 



3 TESTS 

We have tested the algorithm by numerically integrating or- 
bits in a realistic Galaxy potential and after each time-step 
using the above algorithm to determine (8 r ,8 z , J r ,J z ). Any 
variation in the recovered values of the actions along the or- 
bit quantifies errors in the procedure, as do deviations of the 
motion in the (9 r ,9 z ) lane from straight lines. The adopted 
potential is that of model 2 of Dehnen & Binney (1998) 
modified to give the thin disc a scale height of 0.3 kpc - this 
potential is generated by exponential thin and thick stellar 
discs, plus a gas disc, an axisymmetric bulge with axis ratio 
0.6 and a dark halo with axis ratio 0.8. The upper panel of 




0.4 

0.3 0.35 0.4 0.45 0.5 
J, 




0.2 0.4 0.6 0.8 

fl s /27T 



Figure 1. Top: values of J r and J z recovered along an orbit 
in a realistic Galactic potential. The black points are obtained 
with the algorithm of Section 2 using A = 3.5 kpc while the red 
points are obtained with the adiabatic approximation. The units 
are 100 km s kpc. Bottom: the evolution of the angle variables 
along this orbit. 

Fig. 1 shows values of the actions along an orbit that has 
corners at (R, z) — (9.5, 2) kpc and (6.6, 1.35) kpc. The black 
points are obtained using the above algorithm, while the red 
points are obtained with the adiabatic approximation in the 
superior formulation of Schonrich & Binney (2012). Quanti- 
tatively, with the adiabatic approximation the standard de- 
viations of J r and J z are (4.13, 3.89) kms -1 kpc while with 
the above algorithm they are (1.16, 0.97) kms -1 kpc, smaller 
by a factor ~ 4. The lower panel shows the values taken by 
(9 r ,9 z ) at each integration step. The points lie on straight 
lines as required and the slopes of plots of 9i versus time 
agree accurately with the frequencies that are recovered from 
the formulae of Section 2. 

Since the upper panel of Fig. 1 shows that the actions 
we recover, either by the present algorithm or from the adi- 
abatic approximation, are tightly correlated, it is natural to 
ask what else they are correlated with. Their correlations 
with R and z prove to be extremely small (especially in 
the case of the present algorithm), but the red squares in 
Fig. 2 show that in the case of the adiabatic approximation 
Jz (and therefore J r also) is correlated with the combina- 
tion of angle variables 29 r — 9 Z . This angular dependence 
implies that as one moves over an orbital torus at constant 
radius, the error in J z has one sign in the plane and another 
far from it, and that the magnitude of this pattern of er- 
rors oscillates between pericentre and apocentre, changing 
sign somewhere in between. The black triangles in Fig. 2 
show that the present algorithm yields more accurate ac- 
tions largely by eliminating this angular dependence. 

Fig. 3 plots the ratios of the standard deviations of J r 
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Figure 2. J z in units of 100kms~ 1 kpc versus the combination 
of angle variables 28 z —8 r along the orbit that gives rise to Fig. 1. 
The black triangles are obtained with the algorithm of Section 2 
while the red squares are obtained with the adiabatic approxima- 
tion. 




Figure 3. The ratios of the standard deviations in J r and J z 
to (J r + Jz)/2 as functions of the maximum distance from the 
plane attained on the orbit. Along this sequence of orbits J z 
rises from zero to 240kms kpc, while J r decreases from 50 to 
25 km s — 1 kpc. 



and J z to (J r + J z )/2 as functions of the maximum height 
Zmax attained on the orbit - all orbits were started by drop- 
ping particles from (R,z) = (9.5 kpc, z max ). The fractional 
error in J z is never more than 4% and is rarely in excess of 
2%. The error in J r is larger but is still generally under 2% 
of the average action. The pronounced peaks in the errors in 
both actions around z max = 2.5 kpc is probably connected 
with the 1 : 1 resonance between the horizontal and ver- 
tical motions: none of the orbits contributing to the figure 
appears to be actually trapped, but for z max ~ 2.6 kpc the 
frequency Q r — Q z is very low. Consequently, the small dif- 
ference between $ and a Stackel potential has appreciable 
time to disturb the orbit. 

The results shown in Figs. 1 to 3 were obtained with 
A = 3.5 kpc. Fig. 4 shows the standard deviations of J r 
and J z along two orbits as functions of A. The orbits have 
similar eccentricities, but different values of z max : the upper 
squares and triangles are associated with an orbit that has 
z max = 2 kpc, while the lower triangles and points are for 
an orbit that has 2 max = 1 kpc. Both orbits have corners at 
R ~ 9.5 and ~ 6.5 kpc. We see that the standard deviation 
in the values of J z along the orbit is much less sensitive 





Figure 4. The points above SD/J = 0.03 show standard devia- 
tions of J r (triangles) and J z (squares), normalised by ( J r + J z )/2, 
along the orbit that yielded Fig. 1 as functions of the value 
of A used in the algorithm. The lower points show the cor- 
responding numbers for an orbit that has its outer corner at 
(i?, z) = (9.5, 1) kpc rather than (9.5, 2) kpc. The latter orbit has 
actions (45.7, 15.6) kms" 1 kpc. 



to the value of A than is the standard deviation of the J r 
values. 



4 CONCLUSIONS 

We have shown that values of actions and angles accurate to 
a couple of percent can be obtained for orbits in a realistic 
axisymmetric model of the Galactic potential by treating the 
potential as if it were a Stackel potential. For orbits typical 
of observed stars belonging to either the thin or thick discs 
the error in J z is always less than ~ 4% of the average 
action and is usually significantly smaller. The errors in J r 
are always less than 6% and usually less than 2% of the 
average action. Even in the era of Gaia it is unlikely that 
the errors in the measured phase-space coordinates of any 
star will be small enough that the inaccuracies inherent in 
our algorithm will dominate the final uncertainties in derived 
angles and actions. The errors in actions obtained from the 
adiabatic approximation are larger by a factor ~ 4 for thin- 
disc stars and significantly larger still for thick-disc stars. 

A possibility that we have not pursued, but which might 
be important if one needs to model an entire galaxy rather 
than the extended solar neighbourhood, is to make the inter- 
focal semi-distance A a function of L z and E - by integrating 
a few orbits at wide-ranging values of L z and E it should be 
possible to choose a suitable functional form for A(L Z , E). 

Each action evaluation requires a one-dimensional inte- 
gral and with the existing code takes ~ 100 fis on a laptop. 
Each angle evaluation takes about twice as long because it 
requires of order two one-dimensional integrals. Since evalu- 
ation of the observables that follow from a df requires a 
great many evaluations of the actions, it is cost-effective 
to tabulate (J r ,J z ) as functions of the classical integrals 
(L z ,E,Iz) and we have described an effective scheme for 
doing this. In a companion paper we illustrate what can be 
achieved using this scheme by fitting dfs to observational 
data for our Galaxy. 



© 2012 RAS, MNRAS 000, 1-5 



Actions for axisymmetric potentials 5 

REFERENCES 

Binncy J., 2010, MNRAS, 401, 2318 (BIO) 

Binney, J., 2012, arXivl202.3403 

Binney J., McMillan P.J., 2011, MNRAS, 413, 1889 

Binney J., Tremaine S., 2008, "Galactic Dynamics", 

Princeton University Press, Princeton 
Dehnen W., Binney J., 1998, MNRAS, 298, 387 
De Zeeuw P.T., 1985, MNRAS, 216, 273 
Henon M., Heiles C, 1964, AJ, 69, 73 
Kaasalainen M., Binney J., 1994, MNRAS, 268, 1033 
McMillan P.J., Binney, J., 2012, MNRAS, 419, 2251 
Ollongren A., 1965, ARA&A, 3, 113 
Rowley G., 1988, ApJ, 331, 124 
Schonrich R., Binney J., 2012, MNRAS, 419, 1546 



© 2012 RAS, MNRAS 000, 1-5 



